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Abstract 

O 

$h ' We present a detailed numerical study of the equilibrium and non-equilibrium 

dynamics of the phase transition in the finite-temperature Abelian Higgs 
model. Our simulations use classical equations of motion both with and 
, without hard-thermal-loop corrections, which take into account the leading 

quantum effects. From the equilibrium real-time correlators, we determine 
the Landau damping rate, the plasmon frequency and the plasmon damping 
rate. We also find that, close to the phase transition, the static magnetic field 
correlator shows power-law magnetic screening at long distances. The infor- 
mation about the damping rates allows us to derive a quantitative prediction 
for the number density of topological defects formed in a phase transition. 
We test this prediction in a non-equilibrium simulation and show that the 
relevant time scale for defect formation is given by the Landau damping rate. 
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I. INTRODUCTION 

Whilst there are many useful techniques for studying the equilibrium properties of finite- 
temperature field theories, understanding the non-equilibrium dynamics is a much harder 
task. Nevertheless, it would be essential for many fields of physics, for instance cosmology, 
heavy ion physics and condensed matter physics. In all these fields, new empirical data will 
be available in near future, which would allow the theories to be tested, but the complexity 
and the non-equilibrium nature of the phenomena make it difficult to derive theoretical 
predictions that could be compared with the data. 

One fairly generic consequence of phase transitions is formation of topological de- 
fects 0,3. If the phase transition is associated with a spontaneous breakdown of a global 
symmetry, this process is well understood. The correlation length of the order parameter 
cannot keep up with its equilibrium value, which diverges at the transition point. The di- 
rection of the symmetry breaking must therefore be uncorrelated at long distances, and at 
places where these correlated domains meet, topological defects are formed. This is called 
the Kibble-Zurek (KZ) mechanism (see e.g. for a review). 

If the symmetry that gets broken is a local gauge invariance, the above argument cannot 
be used directly, because the direction of the order parameter is not a gauge invariant 
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quantity. We studied this recently in the context of the Abelian Higgs model |[|, and 
pointed out that the thermal fluctuations of the magnetic field lead to another mechanism 
that forms topological defects. The argument was based on fairly generic assumptions, but 
leads to some concrete predictions that were confirmed in numerical simulations. 

The aim of this paper is to study in more detail the dynamics of the Abelian Higgs model 
during the phase transition from the Coulomb phase to the Higgs phase. In particular, we 
concentrate on those degrees of freedom that are relevant for defect formation. This allows 
us to test the scenario of Ref. |4| on a more quantitative level. 

The theory considered in Ref. || was classical, and although the same arguments apply 
to the quantum theory as well, the details of the dynamics are different. The full quantum 
field theory cannot be simulated in practice, but one can argue that the dynamics of the 
relevant long- wavelength degrees of freedom are classical ||. By integrating out the short- 
wavelength fluctuations perturbatively, one obtains a classical effective theory with non-local 
interactions to which we refer as the hard-thermal-loop (HTL) improved theory. In 

order to understand how the quantum effects change the dynamics, we simulate this HTL 
improved theory using the method developed in Ref. ||10|| . 

The structure of the paper is the following: in Sect. |H]we present both the classical and 
HTL improved Abelian Higgs models. In Sect. |TJ we discuss defect formation in the model, 
comparing the mechanism presented in Q with the Kibble-Zurek scenario. In Sects. ^ 
and [V| we describe our numerical simulations and present the results. Conclusions are given 
in Sect. [VI| and technical details of the HTL improved equations of motion and the lattice 
formulation in the two Appendices. 



II. ABELIAN HIGGS MODEL 

The Abelian Higgs model is defined by the Lagrangian 

C = ~F, V F^ + |ZVf - m 2 |0| 2 - A|0| 4 , (1) 

where = <9 M + ieA^ and F^ = d^A u - d u A^. 

A particularly interesting feature of this theory is the existence of Nielsen-Olesen vortex 



solutions |TTfl . These string-like topological defects are characterized by a zero of the Higgs 
field at the centre of the vortex around which the Higgs phase angle has a non-zero winding 
number 



n c 



f dr- V7(f) ^ 0. (2) 
J c 



Here C is a closed path around the vortex and 7 is the Higgs phase angle, i.e., = |0| exp(i7). 



At finite temperature, perturbation theory is plagued by infrared divergences [ I2fl , which 
can be partly cured by a resummation of the perturbative expansion, but even the resummed 
expansion breaks down near the transition. Static equilibrium quantities, such as the phase 
diagram of the theory, can still be calculated reliably with non-perturbative Monte Carlo 
simulations. 

The model has a phase transition between the high-temperature Coulomb phase and the 
low-temperature Higgs phase at T 2 = T 2 « 12(— m 2 )/(3e 2 + 4A). In the perturbative regime 
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(A <C e 2 ), the transition is of first order JDJ, and at larger A it becomes continuous |T4 . 



There are no local order parameters, but a number of non-local ones: the photon mass 
and the vortex tension are non-zero in the broken phase and vanish at the transition |T5 



Therefore the transition is not a smooth crossover like the electroweak phase transition |T6| . 

Monte Carlo simulations cannot be used for real-time quantities in the quantum theory, 
because the necessary path integral is not Euclidean but consists of a complicated path in 



complex time |T7[. However, we can utilize the fact that modes with different momenta 
behave in very different ways ||. The soft, long- wavelength modes (k <C T) have large 
occupation numbers, and they can be approximated very well by a classical theory. This 
makes numerical simulations feasible, because the time-evolution of a classical field theory 
can be found simply by solving the equations of motion numerically. 

A. Classical theory at finite temperature 

Classical field theory at finite temperature is ultraviolet divergent, and thus the results 
depend on the lattice spacing Sx. Divergences like these are generic to all low-energy effective 
theories, and are exactly cancelled by corrections the high-momentum modes induce to 
the effective Lagrangian. If one is only interested in static equilibrium quantities, these 
corrections can be calculated in the limit of high temperature and small lattice spacing 

5x mm, 



2 2 ,2 .^(T 2 3.176T\ 
"* = *" +< 3e +4A »(l2-wj- < 3) 

The approach we use in Sect. [TV] is to take this correction into account and solve the classical 
equations of motion 

dpF 1 ™ = -2eIm0*D>, 
D^D^ = -m 2 T (j) - 2A(>*0)0- (4) 



However, as was pointed out in Ref. [20], these equations do not reproduce the real-time 
dynamics of the quantum theory correctly. Thus, the results are not reliable on a quantitative 
level, but they can still give a reasonably good qualitative picture of the dynamics, and 
provide a non-trivial test for the scenario presented in Section p| . 

B. HTL improved theory 

If the couplings are small, the system is close to thermal equilibrium and T5x ^> 1, it 
is possible to construct a classical theory which approximates the dynamics of the original 
quantum theory to leading-order accuracy in the coupling constants [|]. 

Near the phase transition, m 2 ~ — e 2 T 2 , and we can use the high temperature approxima- 
tion T ^> m in our loop integrals, provided e is small. We calculate the one-loop corrections 
from the hard modes to the self-energies of and A iy and resum them into the effective 



Lagrangian [21 
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+|D„0| 2 - m||0| 2 - AM 4 , (5) 

where is given by Eq. (|3]), and the integration is taken over the unit sphere of velocities 
v = (1, v), v 2 = 1. The Debye mass has the value m 2 D = |e 2 T 2 + 5mfj, where 5m 2 D is a 
counterterm that cancels the UV divergences and is discussed in more detail in Appendix [B 2| . 

All the degrees of freedom remaining in this effective theory are classical, and therefore it 
can be treated as a classical theory. The time evolution of the fields will then be determined 
by the equations of motion [cf. Eq. (§)]: 

9 ,F^ = mlf - 2ehn<f>*D»<f>, 

J Air v ■ a 

D^D^cf) = -m 2 T (p - 2A(0*0)0. (6) 

As such, the equations of motion @ are not well suited for our purposes. Firstly, it is not 
obvious how to find the corresponding Hamiltonian, which is necessary for preparing the 
initial configurations. And secondly, the equations of motion are non-local both in space 
and time. 

In Ref. [fT0|1 , a convenient local formulation was presented along the lines of Ref. ||. 



The latter work introduces a new local field W(t,x,v), representing the departure from the 
equilibrium distribution function for hard particles with velocity v. Since the hard particles 
move at the speed of light, v consists of two free coordinates, making W a six- dimensional 
field. In our formulation, part of the velocity dependence decouples from the dynamics of 
the soft modes, and we can describe exactly the same dynamics with two five- dimensional 
fields, f(t,x,z) and 9(t,x,z), where z G [0,1] is the cosine of the angle between the hard 
mode velocity and the gradient of the distribution function. The details of this formulation 
are given in Appendix |A]. 

In order to simulate the model numerically, we define the theory on a periodic spatial 
lattice, and the z dependence of / and 9 is discretized by expressing them as sums over a 



finite number of Legendre modes, whose order ranges from to iV max ||10|| . The resulting 
equations of motion are detailed in Appendix [A]. We merely note here that, as shown in 
Ref. [ID|, for a given value of iV max and a given momentum k, the approximation breaks 
down at times 

t>t c (k) = W max /k. (7) 

Therefore one can strictly speaking only measure correlators up to At ~ iV max 5x, but since 
the modes that are in equilibrium will remain in equilibrium even beyond that, one can 
essentially trust the results until At ~ t c (k), where k is the relevant momentum scale. 

In Sect. [V[ we study the dynamics of the model using the HTL improved equations of 
motion, and describe how these corrections change the results from the classical case. 



III. DEFECT FORMATION 

In cosmology, temperature of the universe decreases as a result of its expansion, and this 
leads to phase transitions. If we write the equations of motion in conformal coordinates, the 
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effect of the expansion can be absorbed completely into a varying zero temperature Higgs 
mass term, m 2 = m 2 (t). 

In order to see this, we perform the conformal rescaling dx^ — > adx^, — > av^, — > 
a^A^ and — > a~ l <f), where a = a(t) is the scale factor of the Universe. If we can neglect 
the expansion rate in comparison to the microscopic timescales of the theory, which are 
given by m^ 1 and r%j , the effect of this rescaling on the action is to replace the mass terms 
and m 2 D by (mya) 2 and (mzjo) 2 respectively. The hard modes are assumed to be ultra- 
relativistic, so they will stay close to a thermal distribution with temperature T = Ta~ l . 
Hence (m^a) 2 stays constant, as do the rescaled thermal corrections to the Higgs mass, 
leaving the only time dependence in the parameter m 2 (t) = m 2 a 2 (t). In the following, we 
will assume that the phase transition is triggered in this way, by a mass term decreasing 
below a critical value, but we believe that the qualitative features would not be very different 
if some of the other parameters were changing at the same time. 

When the system enters the Higgs phase, Nielsen-Olesen vortices are formed. In the 
limit e — > 0, this can be understood in terms of the Kibble mechanism When m 2 
approaches its critical value, the Higgs correlation length £ grows, and if the system could 
remain in equilibrium it would eventually diverge at the transition point. However, £ cannot 
grow arbitrarily fast, because at the very least it is constrained by the finite speed of light. 
Therefore it remains finite if the transition takes place in finite time. At the transition point, 
the system consists of correlated domains of size £ determined by the maximum correlation 
length reached. In each of these domains, the phase angle of the Higgs field is chosen 
independently of all others, and this gives rise to frustrations, vortices, where these domains 
meet. Up to a numerical factor, the number of vortices piercing a unit area is 

n = N/A « £~ 2 . (8) 

In practice, the maximal rate of change of the correlation length may be well below 
the speed of light, and this argument can indeed be made more precise by considering the 
dynamics of the system in more detail [ 

In Ref. [|J, we argued that this picture is not complete if e > 0. After all, the phase 
angle of the Higgs field is not gauge invariant, and therefore arguments based on it cannot 
apply. Nevertheless, if the amplitude of the magnetic field is small, we can fix the gauge 
in which A4 ~ 0, and then we can use the above picture in this gauge. In these cases 
the above Kibble- Zurek scenario should work. However, if the initial state is at a non-zero 
temperature, magnetic field is never exactly zero, and it must be taken into account. 

In the symmetric phase, the thermodynamics of the gauge field is described by ordinary 
electrodynamics. The energy of any magnetic field configuration is approximately given by 

H EM [B(x)] = ^Jd 3 xB(x) 2 , (9) 

and the probability with which the thermal fluctuations can generate that configuration is 
proportional to exp(—HEAi/T). Let us consider a circular loop C that bounds a surface S. 
The magnetic flux through the surface is 

$ 5 = J dS-B. (10) 
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Although this is zero on the average, it is non-zero in almost all configurations, and therefore 
the typical value, given by is non-zero. We can estimate its value by calculating the 
energy of the field configuration that minimizes the energy for a given value of This 
minimal configuration is simply a magnetic dipole, and its energy is 

£min($s) ~ (11) 

where R is the radius of the loop C. Thermal fluctuations can create this configuration if 
T> E min (<& s ), and solving this for $5 shows that the typical flux through the loop is 

$5 « VTR. (12) 

When the system enters the Higgs phase, magnetic flux is confined into flux tubes, which 
costs energy. Therefore the dynamics tries to decrease the magnetic flux, but it can only do 
so at short distances, because the range of the interaction between the flux tubes decreases 
rapidly. 

We can be more specific in Fourier space. We can write the two-point correlator of the 

— * 

magnetic flux density B as 

(B^B^k')) = ^, - ^ (2nf5(k + k')G(k). (13) 

In the symmetric phase, different Fourier modes B(k) behave as independent oscillators in 
thermal bath, and thus each of them has the same amplitude Go(k) = T. When the system 
enters the broken phase, magnetic field becomes massive, and the equilibrium distribution 
changes into 

k 2 

G «- T FT< (14) 

In order to stay in equilibrium, the amplitude of the long- wavelength modes must drop 
rapidly, but the time scale r(k) of the dynamics of these modes is very slow. It depends on 
the individual system, but in general limfc^ r (^) = 00. Thus the modes with k less than 
some critical value k cannot remain in equilibrium. If we know r(k), we can calculate k 
from the condition 



d\nG(k) 



dt 



r(k) 



(15) 



The consequence of the above process is that there will be long-wavelength magnetic 
fields present even in the Higgs phase. In particular, at distances less than £ = 2n/k it looks 
like there was a uniform external magnetic field. We can estimate that its amplitude _B avg is 

w? Bfk) r Sr B0) ) x r S g ° w ~ rp - <i6) 

This magnetic flux must be confined into vortices, and because each vortex carries one flux 
quantum $ = 27r/e, the number density of vortices per unit area is 
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(17) 



Even if we do not know r(k), we can still make some concrete predictions based on 
this scenario. For instance, the spatial distribution of vortices turns out to be completely 
different than in the KZ scenario. At distances shorter than £ = 2%/k, the magnetic field 
points in the same direction, which means that the vortices tend to be aligned, while in the 
KZ scenario they prefer to be anti-aligned. This prediction was confirmed in the simulations 
in Ref. §. 



IV. CLASSICAL SIMULATIONS 
A. Equilibrium 



In order to test the scenario presented in Sect. fTJ, we carried out a number of numerical 
simulations. Let us first discuss the simulations of the classical model (§). 

Because we are interested in transitions that start from close to thermal equilibrium, 
we will first have to thermalize the system. This means preparing an ensemble of field 
configurations with the canonical equilibrium distribution exp(—/3H), where (3 = 1/T and 
the Hamiltonian H is 



H= d 3 x 



]-E 2 + UV X Af + 71*71 + (A0)*(A0) + "40*0 + A(0*0) 2 



(18) 



Here tt = <9o0 is the canonical momentum of 0, and in the temporal gauge (Aq = 0), the 
electric field is simply E = — 8qA. In addition, the fields must satisfy the Gauss law as an 
extra constraint 

V • E = 2elm0*7r. (19) 

Because of the constraint (|T9|), a straightforward Metropolis algorithm would not work 
very well. Instead, we used a hybrid algorithm, in which we thermalized the component of tt 
orthogonal to Eq. flT9"|) with a heat bath algorithm, and performed a number of Metropolis 
updates to the gauge field A. Because A does not appear in Eq. (|19"D, it leaves the constraint 
unchanged. Then we evolved the system with the equations of motion 

d A = -E, 
<9 O = vr, 

d E = V x V x A + 2elm0* A0, 

d 7r = DiD4 - m|0 - 2A(0*0)0, (20) 

which, again, leave the Gauss law unchanged. We repeated this procedure a number of times 
so that the system thermalized. 

In our simulations, we used the couplings e = 0.3 and A = 0.18. The lattice size was 
120 x 120 x 20 (the reason for choosing one short dimension is discussed in Sec. [IV B| ), the 
lattice spacing was 5x = 6T _1 and the time step was 5t = 0.3T -1 . The details of the lattice 



implementation are given in Appendix B 1 
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Since A > e 2 , the phase transition is continuous. In order to determine the location of 



the transition point, and to test the accuracy of the tree-level result (14), we measured the 
correlator ( |T3| ) at various values of m 2 , starting deep in broken phase. For thermalization to 
each value of m 2 , we used 24 hybrid Monte Carlo cycles each consisting of 400 Metropolis 
sweeps and time evolution for At = 600T -1 . We carried out the measurement in nine inde- 
pendently thermalized configurations, measuring the average correlator in an unperturbed 
run of length At = 12000T" 1 (and At = 72000T- 1 for m 2 = -0.083T 2 ). The results are 
shown in Fig. |l]a. The solid lines show that the agreement with Eq. ([14] ) is excellent, and 
we find that the phase transition takes place at m 2 « — 0.083T 2 . 

In the symmetric phase, the tree-level result fll4]) corresponds to a constant G(k) = T, but 
the data measured at m 2 = — 0.083T 2 (see Fig. ||b), clearly turns down at small momenta. 
It is customary to parameterize the corrections to the tree-level result by introducing the 
static photon self energy U T , defined by 

Our measurements seem to contradict the results of Kraemmer et al. and Blaizot et 
al. |21|j2^|, who showed that after a resummation, the lowest-order term of 11^ is proportional 



to k 2 . Such a quadratic term would only change the overall normalization of the correlator 
and any higher-order terms would only modify the high-fc end of the spectrum. However, 
they considered specifically the case in which the zero-temperature Higgs mass vanishes, i.e., 
m 2 = 0. In that case, the thermally generated effective mass for the Higgs field M 2 , which is 
approximately equal to m 2 -, is always ~ e 2 T 2 [see Eq. (^[)] and acts as an infrared regulator 
in the loop integral, making Hr(/c) an analytic function of k 2 . Because one can show that 
the constant term is forbidden, the lowest-order term must be 0(k 2 ). 

In our case, this argument breaks down because we are studying the system at the 
transition point where the effective Higgs mass M becomes very small. In the momentum 
range M <C \k\ <C eT, the perturbative behaviour of Hr{k) at one loop order is 

e 2 T 

U T (k) = —k = 0.005625T&. (22) 

This is shown as a solid line in Fig. |l|b, and agrees fairly well with the data. Assuming 
that perturbation theory is applicable, the discrepancy at very low k can be explained by 
a non-zero effective Higgs mass, which arises because we are not exactly at the transition 



point. With M > 0, the one- loop result is |£3 

e 2 MT 



n r (fc) 



47T 



4M 2 + k 2 k 

— — — arctan — — — 1 

2kM 2M 



(23) 



The best fit gives M = (0.0040 ±0.0005)T, and is shown in Fig. |I|b as a short-dashed line. It 
is curious how well it agrees with the measurements, because one would expect perturbation 
theory to break down at low momenta near the phase transition. 

Assuming that the Higgs mass vanishes at the transition point, all that remains is 
Eq. fl2"2D - The existence of a linear term in the self energy Hp is surprising, because it 
implies magnetic screening. This screening is not as strong as it would be if there was a 
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FIG. 1. The spatial correlator of the magnetic field in the Fourier space, a) The solid lines 
are fits to Eq. ( |l4|) and show that the system is in the Higgs phase when m 2 < — 0.083T 2 . The gray 
squares show the spatial correlator after a quench as discussed in Section IV B. The open circles 
correspond to the HTL simulations in Section V A , b) The m? = — 0.083T 2 data on a linear scale. 
The short-dashed and long-dashed lines show the fits to Eqs. ( |23"| ) and ([?(]), respectively, and the 
solid line shows the perturbative result (^). 



constant term, in which case the correlator would fall exponentially in the coordinate space. 
With a linear term, the low-momentum behaviour of the transverse gauge field correlator 
is ~ and consequently, the long-distance behaviour in the coordinate space is ~ r~ 2 
instead of the usual ~ r _1 . 

We also measured real-time correlators in the same simulations. In each of the nine 
independent configurations, we measured the correlator G(t, k) and took the average of the 
results. Two examples, measured for k = 0.026T at m 2 = — 0.083T 2 and m 2 = — 0.086T 2 , 
are shown in Fig. ^|a. At all values of m 2 and k we measured, the data were well described 
by the function 

G fi t(t) = a exp(-7 L t) + a x exp(-7 P t) cos(u p t + 5), (24) 

where a , a±, jl, j p , u p and 5 are free parameters. Physically, 7^ is the Landau damping 
rate, 7 P is the plasmon damping rate, and uj p is the plasmon frequency. We fitted this 
function to the data at each point and estimated the errors using the bootstrap method. 
The results are shown in Figs. 0b d. 

We can compare these results with perturbative calculations, but we must keep in mind 
that since these simulations were carried out in a classical lattice theory, the result is not the 
same as in quantum theory in continuum. In the symmetric phase, the plasmon frequency 
should behave as 

oo p (k) = ^ 2 + m 2 , (25) 

where m p is the plasmon mass, which has been calculated in classical lattice perturbation 
theory in Refs. HH, 
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ml « 0.086e 2 -^- « 0.0013T 2 . (26) 

We have plotted this curve in Fig. ^|b, and it agrees very nicely with the measured frequencies 
at m 2 = — 0.083T 2 . In the Higgs phase, the photon becomes massive due to the Higgs 
mechanism and this increases m p . 

The plasmon damping rate has not been calculated perturbatively for the classical lattice 
theory. However, a calculation has been carried out in the HTL-improved theory by Evans 



and Pearson ||24|| , who showed that it was peaked just below the phase transition. Our 
classical results also show a rising trend as the transition is approached, although we cannot 
directly compare the numerical values. 

The Landau damping rate in the symmetric phase has been calculated perturbatively in 
Ref. IE5I, ' 



lL ~ ^ « 470^. (27) 

However, in our results the dependence on k is milder and we find the best fit with jl w 
15T k above the transition point. Below the transition point, the exponential decay 
rate 7^ becomes very large so that at late times, the correlator simply oscillates around zero 
(see Fig. ||a). In Fig. ||d, this behaviour can be seen in the values of the Landau damping 
rate. 

In the defect formation scenario, the freeze-out occurs when Eq. ([15]) ceases to be satisfied. 
This happens most easily at the transition point, and there the real-time correlators are still 
the same as in the symmetric phase. Therefore, it is the symmetric phase correlators that 
determine the relevant time scale r(k), and Fig. |2] shows clearly that for the long- wavelength 
modes, the lowest time scale is that of Landau damping, r(k) ~ •y L (k)~ 1 . 

Substituting r(k) = l L {kY l to Eq. (|TJ), we find 



1 5m 2 
k 2 t q 



7l « 15T- 1A k 2 -\ (28) 



which implies 



'0.006l\ 1/4,1 



T « 0.29T(Ttq)-°' 24 oc tq - 24 . (29) 



Tr< 



Q 



B. Phase transition 

We simulated the phase transition by thermalizing a number of configurations with 
m 2 = m\ = — 0.044T 2 , and solving numerically the equations of motion (20) , varying the 
mass term with time according to 

/ 4 1\ 
m 2 (t) = ml — 5m 2 ( — arctan(t/rQ — 1) H — I , (30) 
\3n " 3/ 
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FIG. 2. The real-time correlators measured at different values of m 2 . a) Two examples 
of real-time correlators, measured at m 2 = — 0.083T 2 (solid) and m 2 = — 0.086T 2 (dashed), 
k = 0.026T, together with fits of the form (24). b) The fitted plasmon frequencies u p , and the 



perturbative estimates, c) The fitted plasmon damping rates j p . d) The fitted Landau damping 
rates 7^, and the perturbative estimates. In plots b-d, the gray squares correspond to the state of 



the system after a quench and open circles to the HTL simulations discussed in Section V A 
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where 5m 2 = 0.089T 2 . This form of m 2 (t) has the advantage that long time after the 
transition, the system is in thermal equilibrium, making it possible to compare the final 
states. 

The vortices produced in the transition are closed loops, and in general they will soon 
shrink to a point and disappear, which makes it very difficult to even define what we mean by 
the final vortex number. However, on a periodic lattice some of the vortex loops can be non- 
contractible, i.e. wind around the lattice in some direction. These loops can only disappear 
if they annihilate with another vortex of opposite direction, but this process is very slow, 
since the interactions between the vortices are exponentially suppressed at long distances. 
Therefore we chose one of the lattice dimensions much shorter than the other two (but still 
much longer than the microscopic length scales such as the Debye screening length m^ 1 .) A 
long time after the transition, when the system is deep in the broken phase, these vortices 
still remain and are well-defined macroscopic objects, and they can be counted without any 
ambiguities using the gauge-invariant lattice definition for the winding number fl2"7| , |2"§f . 

Because the scenario of defect formation discussed here is based on the assumption that 
the distribution of the magnetic field freezes in the transition, we can carry out a very 
simple test for the scenario by quenching the system through the transition and measuring 
the spatial correlator. We used tq = SOOT^ 1 and stopped the quench at m 2 = — 0.089T 2 to 
carry out the measurement. The spatial correlator is shown in Fig. |l]a as gray squares, and 
indeed, resembles much more the symmetric phase correlator than the equilibrium correlator 
at the same value of m 2 . We also measured the real-time correlator, and the corresponding 
time scales are shown in Figs. 0b-c. The plasmon frequency and decay rate do not differ 
significantly from their equilibrium values, but Landau damping gets extremely slow. This 
means that once a mode has fallen out of equilibrium at the transition point, it takes a very 
long time before it thermalizes. 

When one of the dimensions is shorter than £ = 2ir/k, the prediction (|i"7|) changes into 0] 

n « —T^L-^k- 1 . (31) 

We tested this prediction by simulating the time evolution with different values of tq and 
measuring the vortex number a long time after the transition at t = tq + 2400T -1 . The 
results were published in Ref. ||, and are shown in Fig. [| Each point is an average over 
around 15 runs starting from different initial configurations. 

Combining the result k ~ t q°' 24 from Eq. (|29|) with Eq. (|3T|), we expect iV « ctq ' 2 . 
If we determine the constant c from the best fit to the data at tq > 200T -1 , we find 
c = 39.5 ±1.0 and x 2 = 7.6/13 dof. This is shown in Fig. |^ as the dashed line. If we 
also leave the exponent as a free parameter, we find N = (43.9 ± 7.9)(Ttq) _0 ' 255±0 ' 026 , with 
X 2 = 7.2/12 dof. 

If we calculate precisely the prediction ( |3~1~1) using Eq. (p9|) , we find c ~ 650 for the 
above fit parameter. However, it is not surprising that it differs from the measured value 
by a factor of « 15, not only because have we neglected factors or 2tt and other numerical 
factors, but also because in our estimates we have assumed that the mechanism is ideally 
efficient and that all the magnetic flux is converted into vortices. 

From the cosmological point of view, we are more interested in the area density of vortices 
in a fully three-dimensional case than in a thin box. Because the topology of the system 
does not prevent vortex loops from shrinking, the resulting network is not stable, unless 
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100 1000 10000 

FIG. 3. The dependence of the final number of vortices on the quench rate tq. The dashed 
line is a power-law fit of the data at tq > 200T -1 with the exponent —0.24 predicted by Eq. (|29|). 
The circles and squares correspond to the HTL simulations discussed in Section [VB| , and the solid 
line is a power-law fit of tq > 100T -1 using the exponent —0.2 predicted by Eq. ([42]) ■ 



it is stabilized by the expansion of the universe. We can therefore only use Eq. (|T7| ) as 

an estimate for the area density of vortices immediately after the transition. Comparing 
Eqs. flT?P and (0), we find 

%D=(|) 1/2 T^/ 4 ^ 2 , (32) 

where is given by Eq. ([H]) and n^n is the three-dimensional area density given by 
Eq. (13). In our case, L z = 120T" 1 and A = 5.2 x 10 5 T" 2 , and we find 

n 3D « lQQT^nlg « 1.3 x 10-\Tt q )-°- 38 . (33) 



V. HTL SIMULATIONS 
A. Equilibrium 

As discussed in Section O, the classical theory discussed above does not describe the 
dynamics of the quantum theory correctly. Therefore we also carried out the same simu- 



lations with the HTL improved theory fl5|). In the same way as in Section [IV A[ we first 
studied the equilibrium properties of the theory. We used the same couplings e = 0.3 and 
A = 0.18 as in the classical case, and the same lattice. With these parameters, the Debye 
mass has the value m 2 D = 0.03T 2 . The number of Legendre modes was iV max = 4. The 
only effect of the HTL corrections to the thermodynamics is to give an extra contribution to 
electric screening, so we expect that the phase diagrams of the theories are practically the 
same. Therefore we only carried out the equilibrium measurements at the transition point, 
m 2 = -0.083T 2 . 

Using the Hamiltonian (|A"8|), we first prepared the thermal initial conditions with a 
Monte Carlo algorithm. This can be done in two steps: 
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1. Generate the soft mode configuration. The part of the Hamiltonian that involves the 
hard modes is Gaussian and therefore they can be integrated out exactly. This results 
in a theory with only the soft fields and and their canonical momenta. The Gauss 
law leads to an extra Debye screening term 

5H Dehye = J_ . £ - 2elm0*7r) . (34) 

in the Hamiltonian.Q We used a Metropolis algorithm for and a heat bath algorithm 
for all other fields, and carried out around 10000 thermalization sweeps. 

2. Generate the hard modes in the background of the soft modes generated in step 1. 
Since the Hamiltonian is Gaussian it could in principle be diagonalized and the field 
values could be taken directly from a normal distribution. However, we did this only 
for 9 and its momentum n. Note that the lowest Legendre mode of n is fixed by the 
Gauss law, 

n (0) = — (V • E - 2elm0*7r 

m D v 

The field / and its momentum F were generated using a heat bath algorithm with 
around 5000 sweeps. 

Again, we evolved the configurations taken from the thermal ensemble using the equa- 
tions of motion ( |A6| ) for the time At = 12000T -1 , measuring the equal-time and real-time 
correlators. The equal-time correlator is shown as open circles in Fig. [I], and it agrees 
fairly well with the corresponding classical correlator. This time, the fit to the perturbative 
one- loop result (p3|) favours M — 0, which suggests that because of the presence of the 
hard modes, the transition point is at slightly larger m 2 than in the classical case, and that 
m 2 = — 0.083T 2 is very close to it. 

Because perturbation theory cannot be trusted at low momenta, we do not adopt the 
perturbative result ([22]), but instead we simply assume that the self energy is linear at small 
momenta 

n T (fc) = k np k, (36) 
and determine the coefficient k np from the best fit to the data, which gives 

k np = (0.0071 ± 0.0005)T, 
and use this in our later estimates. 



1 In principle, the canonical momenta could also be integrated out at this stage, resulting in 
a theory with an extra neutral scalar field Aq. This would lead to an effective theory that is 
equivalent with dimensional reduction at one-loop order |29| . 
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(35) 



(37) 



The hard modes have a significant effect on the real-time correlators, shown in Figs. |2|b- 
d. Because they mimic the effect of the hard modes in the continuum quantum theory, we 
should now be able to compare the results with the standard perturbative calculations. For 



the plasmon mass, the perturbative result is |30| . [21 

1 



m 



P %/3 



-^m D = 0.1T. (38) 



The dashed line in Fig. |2|b shows the corresponding curve u p = ^Jk 2 + m 2 . The measured 
values are slightly below this curve, but the agreement is still very good. 

The continuum plasmon damping rate at zero momentum has been computed perturba- 



tively in the HTL approximation in Refs. [[21| ,f24 

l P (T) = ^A(T/T c ,X/e 2 ). (39) 

The function A(T/T C , A/e 2 ) takes the value 1 at the critical point, peaks at value ~ 1 just 
below it, and vanishes above it when m p > 2M. In our case, the damping rate at the 
transition point would be 

7 P (T) « 0.0012T, (40) 

which is shown in Fig. |2]c as a dashed line and is lower than value the measured at m 2 = 
— 0.083T 2 by roughly a factor of three. It has been observed earlier ||26|| , that also in the SU(2) 



theory, the perturbative calculation underestimates the plasmon damping rate significantly. 
We can also note that the damping rate is lower than in the classical theory. 

The Landau damping rate can be obtained directly from the HTL improved Lagrangian 



and is 30 



4k 2 k 2 
1l = 2~(k + k np ) « 42.4— (k + k np ), (41) 

7T7TI jj J. 

where we have taken into account the linear contribution ([36]) to the self energy. This curve 
is shown as the dashed line in Fig. |2|d, and agrees well with the measured values. Again, we 
can conclude that the relevant time scale is that of Landau damping. If we can ignore the 
linear correction, we find [cf. Eq. (|29|)] 



V4 u r Q 

and in slowest transitions, we have 



- ( -m 2 D ^— \ ' « 0.29T(Tr Q r - 2 ~ r Q - 2 , (42) 



% ^ ^m|^y 25 ^ 74T(Ttq) -o. 25 ^ r -o. 25 (43) 

However, this corresponds to tq ^> 10 8 T _1 , which is much larger than the values of Tq we 
are able to use in out simulations. 
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B. Phase transition 



As in Sect. [IV B| , we studied the non-equilibrium dynamics of the phase transition by 
starting from thermal configurations at m 2 = — 0.044T 2 and evolving the system with a 
time-dependent mass term (|30|). We used two different values for iV max , 4 and 16. 

In Fig. H], we show how the final vortex number, measured at t = 2tq + 2400T~ 1 depends 
on the quench rate tq. Each data point is an average over 20-30 runs. In fast transitions, 
tq< 1000T _1 , the results for iV max = 4 and 16 agree, but in slower transitions there is a 
statistically significant difference. This can be understood in terms of Eq. (|7]): For small tq, 
-^Vmax = 4 is sufficient because by the time the approximation breaks down, the system is 
already so deep in the Higgs phase that the vortex number cannot change any more. When 
tq gets larger, eventually a point is reached at which the breakdown occurs so early that it 
would have an effect on the final state. This may have happened even for iV max = 16 in the 
slowest transitions with tq = 6000T -1 . 

Combining Eq. (fH with Eq. we find N ~ r Q - 2 . A fit to the iV max = 16 data at 

tq > WOT' 1 with N = c(Tt q )-°- 2 gives c = 45.5 ± 0.9 with x 2 = 6.8/6 dof, which shows 
that the results are compatible with the prediction. The fit is shown in Fig. § as a solid line. 
The prediction of Eq. (|31| ) is c ~ 660, which is greater than the measured value by a factor of 
« 15. This is the same factor found in the classical case, which further supports our scenario. 
If we keep the exponent free parameter, we find N = (45.9 ± 4.6) (Tr Q )- a201±0 - 015 , with 
X 2 = 6.8/5 dof. 

Again, we relate the results to three dimensions using Eq. ([52|), and find 

n m « 1.4 x 10- 4 (t q T)-°- 30 . (44) 

The conjectured non-perturbative behaviour at low momenta ( |3"E| ) would change this scaling 
law in very slow transitions. Eq. (|TTD would become n oc k , and Eq. ( [43]) would therefore 
imply n oc Tq ' 5 . 



VI. CONCLUSIONS 

In this paper we have presented a thorough study of the dynamics of the high temper- 
ature phase transition in the Abelian Higgs model, using both classical and HTL improved 
approximations. The aims were twofold. The first was to measure the equilibrium properties 
of the theory, and in particular find the relevant timescale for the decay of the long wave- 
length modes of the gauge field, which were identified in |4[] as crucial to the understanding 
of vortex formation. 

In Ref. [|J, we performed numerical simulations of the Abelian Higgs model phase transi- 
tion in real time by using the classical theory and changing the mass parameter of the Higgs 
field over a characteristic quench time tq. In this paper, our second aim was to do the same 
quenches with the HTL improved theory, and to compare the resulting scaling law for the 
number of vortices N formed as a function of the quench time. 

Our measurements of the equilibrium correlators show that perturbation theory gives a 
reasonable account of their behaviour, except perhaps for the plasmon damping rate. A new 
and unexpected result is that near the transition, the equal-time correlator exhibits power- 
law magnetic screening, with a coefficient that is similar but not equal to the perturbative 
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one. In the HTL improved simulations, the numerical values of the plasmon mass and the 
Landau damping rate agree well with the perturbative values in the Coulomb phase. For 
the plasmon damping rate, however, the measured value was significantly higher than the 
perturbative estimate. As might be expected, the agreement is not as good in the classical 
theory: although the plasmon mass agrees well, the dependence of the Landau damping rate 
7l on wavenumber k is ~ k 2 ' 1 rather than the expected k 3 . 

We found significant differences between the scaling laws for the classical simulations, 
and for HTL improved quenches with iV max = 4 and A max = 16. The difference between 
the classical scaling law and the HTL improved one with A max = 16 can be ascribed to the 
discrepancy in the Landau damping rate, and lends force to the contention made in that 
it is the balancing of the cooling rate with the Landau damping rate which decides the length 
scale above which the fields fall out of equilibrium. The difference between iV max = 4 and 
A max = 16 can be understood as stemming from a lack of phase space in the hard modes, 
which causes the HTL approximation to break down at times greater than t c (k) = 4iV max /fc. 

Our simulations have been carried out to leading order in the couplings e and A, and 
hence do not include the effect of high momentum transfer scattering between the hard 
modes. However, the hard modes can still scatter by exchanging soft modes, we do not 
expect this to change the dynamics qualitatively on the relatively short time scales which 
we have been able to study. In very slow transitions, it may be that hard mode scattering 
and also non-perturbative effects in the photon self-energy start to become important and 
and the predicted scaling law n oc Tq ' 3 ceases to be valid. 
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facilities, funded by HEFCE, PPARC and SGI. We also acknowledge computing support 
from the Sussex High Performance Computing Initiative. AR was supported by PPARC 
and also partly by the University of Helsinki. 

APPENDIX A: HTL IMPROVED EQUATIONS OF MOTION 

In Sect. the local formulation of the HTL improved Abelian Higgs model was described. 
In this Appendix we detail the resulting equations of motion for the fields <fi and A and 
Legendre modes 6 and /, which encode the effect of high momentum (k > T) particles. 

These fields satisfy the equations of motion 
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-V x V x A-2ehn<j)*D<j) 




(Al) 
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Since these equations for / and 9 are linear, it is easy to solve them and show that the 
dynamics of <fi and Ai is identical to the original non-local theory (H). 

Not only is this reformulation of the theory local, but the equations of motion are in a 
canonical form and we can therefore write down the corresponding Hamiltonian 



H 



d x I dz 
o 



1^2 1 



-E z + -(V x A) + ir*n + (A0)*(A0) + m l T (j)*(j) + \{<t> 



2 2 

,2 



+\f 2 + hi 2 + Z -{V x f) 2 + Z -{Vd - m D Af - m DZ] 



1-z 2 



-f-VxA 



(A2) 



where F = dof and II = 8q9 are the canonical momenta of / and 9, respectively. We also 
need two extra conditions, namely the transverseness of / and Gauss's law 



V-/ = V-F = 0, 



V • E = -m D / dzU(z) + 2elm0*7r 
Jo 

The z dependence of / and 9 is discretized jLOj by introducing the Legendre modes 

fin) 



(A3) 
(A4) 



1 dzzJ-^—P 2n (z)f(z), 0< n > = f 1 dzP 2n (z)9(z), 
o V 1 — z Jo 



Fin) = / ^2n(z)F(z), nW = / dzP 2n (z)U( Z ). 

Jo z V 1 — z Jo 



(A5) 



Note that we have used slightly different definitions from Ref. ||10|| , in order to write the 
Hamiltonian in a practical form. 

In terms of these modes, the equations of motion become 

d A = -E, 

d 9 (n) = n (n) , 

9 £ = VxVx A + 2elm(p* D(p + -m 2 D A 

3 

~m D (W (0) + 2W (1) + V x /t°) - V x /tD) , 
do pH = _y x y x f\n) + mo f x AS nQ} 

doU (n) = C n + V 2 e (n+1) + C°V 2 # (n) + C; t V 2 9^ - m D V ■ A (^5 n>0 + ^5 n ^ , (A6) 



where 

= (2n + l)(2n + 2) 
B (4n+l)(4n + 3)' 



'(2n + l) 2 ^ 4n 2 



An + 1 \ 4n + 3 4n - 1 



2n(2n - 1) 
(4n + l)(4n- 1)' 

(A7) 



The Hamiltonian becomes 
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1* 1 



H = / d 6 x\ -E + -(V x A) + n*n + (A0)*(A0) + + A(0 



+ E 

n=0 



. 2 2 

co rl 1 



4 4n + 3 



((2n + l)F (n) + (2n + 2)F (n+1 - 



+ 



4n + 1 



n) 



2 1 1 



4 4n + 3 



(2n + 1)V x + (2n + 2)Vx/ 



n+l) 



4n + 1 /„,^\2 1 1 



2 



( n (n) ) + f (2n + 1) W (n) + (2n + 2) V0 (n+i; 

v / 2 4n + 3 ^ 

V • A (9^ + 20«) - x A (/t°) - /«) + im^i 2 



(A8) 



APPENDIX B: LATTICE DISCRETIZATION 
1. Classical theory 

In order to carry out numerical simulations described in Section we discretize the 
Hamiltonian (18) and the equations of motion (|20| ) in the standard leap-frog fashion. The 
Higgs field was defined at the lattice sites, and its canonical momentum at temporal links 
between time slices. The gauge field A was represented by real numbers defined at links 
between lattice sites, and the electric field E by temporal plaquettes that connect these 
links. Therefore 717+^ is actually defined at the point (t + 8t/2,x), Aiug\ at (t, x + i/2) and 
Ei,H,£) a t (t + 8t/2, x + i/2). Here % is a vector of length 8x in the % direction. 



The lattice version of the Hamiltonian ( 18|) is 



H = Y J 8x z 

X 

2 



I E E i + \ E (%feA+^) 2 + 7T*7T 



6 



^Re^f/^)^) + (to 2 . + 



where 



:bi) 



Ui = exp (ze<5xAj) , 
Af 0(5) = ifoT 1 (0(5±») - 0(5) 

We also define the lattice version of the covariant derivative 

Dt<P(x) = (^i,(5)0(5+i) - </>(*) ) 7 



A 4>(x) = 8x 1 



f.F) 



u, 



(B2) 



(B3) 



The value of the bare lattice mass m 2 - is given by Eq. (|3]) and was chosen in such a way 
that the Hamiltonian ( [Bl~|) describes the thermodynamics of the finite-temperature theory 
with renormalized mass m 2 correctly 
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The discretized equations of motion are 



\Ei,(t,z) — tijktkimAj Ai t (t t g) — 2elm(j)* t ^D i (j)(t,s), 



A<f>(t+St,2) = K(t,x), (B4) 

where A t (t) = <ft -1 [0( t ) ~ <P(t-st)] etc. 

The lattice version of the Gauss law ([19]) is 

E A i E h(t,s) = Selm^^TT^). (B5) 

i 

This is an extra constraint the initial field configuration must satisfy. 



2. HTL theory 



In the HTL simulations described in Section |V|, the soft modes were discretized in the 
same way as in the classical case. The extra field 9 is defined at lattice sites and fi at the 
plaquettes. We denote by fi t (t,x) the field value at (t, x + x/2 + y/2 + z/2 — i/2). 

The lattice version of the HTL-improved Hamiltonian ( |A8[ ) is 

h htl = H + 5x 3 J2[^F + n f + n n + n e ], (B6) 



where 

n 



n=0 ^ 



4n + 3 



((2n + + (2n + 2)F (n+r 



-(4n + 1) (C+F (n+1) + C°F (ri) + C~F {n - r > 



4n + 1 



^n = E 



n=0 

+(2n + 2)ej ifc A7/ fc 
^ 4n + 1 /_,^n2 



_ ,(n+l) x2 ' 



4 4n + 3 

m ^ A+4 fW"l _i_ X m 2 4 

—e ijk /\j A k \f k - f k j + -m D A, 



n=0 
oo 



(n(«))' 



11 9 m 

W fl = £ ((2n + l)A+fW + (2n + 2)A+^)) + ^Ar^ (0<°> + 26lW) . (B7) 



The Gauss law (IA4D can be written in the form 



n 



(0) 
(t,x) 



— ( E A i E i,m - 2eW^7r (t ,^ 



(B8) 



and therefore we can eliminate II^ ^ from the Hamiltonian. 
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The bare Higgs mass has the value given in Eq. (W), and the Debye mass is [18.19 



mfl = 3 eT - 2e W (B9) 

These counterterms were calculated by matching static correlators and in the absence of 
Lorenz invariance, they do not remove all the ultraviolet divergences from real-time quanti- 
ties. However, since our lattice spacing is relatively large, this leads only to small errors. 
The discretized equations of motion are 

\Em,£) = £ijk£kimAjA l Ai t (t,g) — 2elm<j)*( t ^D + <j)i^t,x) 

g - {^i U (t,x) + llX i V(t,x) + tijk^j Jk,(t,x) ~ tijk^j Jk,(t,x) 

MS) = Wm 1 ' + <WM5) + C-AtA~9^ 

-m D Aj Ai^ t ,3) (3^,0 + > 
A t (f)(t+8t,x) = ft(t,x), 

M&mo = n S)- ( B1 °) 
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